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1. Introduction 

Two-point Feynman paiameter integrals with one mass containing local operator insertions in D- 
dimensional Minkowski space with one time- and (D — 1 ) Euclidean space dimensions (e = D — 4 
and £ G M with |£| <C 1) can be transformed symbolically to multi-sums over hypergeometric 
terms depending on e and a discrete Mellin parameter n. A detailed description of the underlying 



algorithms was given in [22]. Sums of this kind emerge in the calulation of Feynman diagrams 



mentioned at 2- and 3-loops [18, 11, 13, jTS], ^ ^] which are of importance for the precision 
measurement of the strong coupling constant and the parton distribuitions from the world deep- 
inelastic data, cf. e.g. [jl^; at the 3-loop level a larger number of moments of these quantities has 
been calculated in [[I|, |0|] ^ 

Given such a sum representation like, e.g. in ^ ^, |l^, |5[ ^ Q], the challenging task is to 
simplify these multi-sums (usually several thousands coming from Feynman integrals that describe 
a physical problem), in terms of well known special functions that can be processed further. In order 



to accomplish this task, a difference field theory for symbolic summation [Eq, 29 , 3 1 , 32 , 33 ] based 



on nZ-fields [ [24| | is heavily used within the summation package S i gma [ )30| ] that generalizes the 
summation paradigms of [^] to multi-sums. In addition, symbolic and analytic techniques for 
harmonic sums [ 34 , [l^, |l^] and generalizations to cyclotomic sums ^ implemented within the 
HarmonicSums package [|]] are exploited. 

In this article we report on new summation technologies that enables us to simplify huge 
multi-sum expressions coming from Feynman integrals to expressions in terms of indefinite nested 
product-sums. This allows us to represent these expressions, if possible, in terms of harmonic 
sums and if this is not possible, in terms of ^-sums or more generally in terms of cyclotomic sums; 
in any case, the result is given in terms of products and sums which are algebraically independent 
among each other [17, ^]. 



The resulting ideas are implemented in the two Mathematica packages EvaluateMultiSums 
and SumP reduction. The procedures are illustrated by a new calculation of fermionic contri- 



butions to the gluonic massive operator matrix elements []23[]. 
2. Evaluation of Multi-Sums 

Given a Feynman integral F{n) depending on a discrete Mellin parameter n and a small dimensional 
parameter e, one is interested in the first s coefficients in its Laurent series expansion 

F{n) = F,{n)e"+F,+ y{n)e"+' + ■ ■ ■ + F,+,^y{n)e"+-'-' + 

M G Z. As worked out in [^], a large class of Feynman integrals can be transformed to a multi-sum 
written in the form 

Li(«) Ly{n,k[,...,ky-i) I 

F{n) = £ ... £ £/(£,«,^i,...,fc„) (2.1) 

kl—ll ky—ly k—[ 

where L,(«,/:i, ... ,^y_i) stands for an integer linear relation in the variables n,ki, . . . ,k^,^i or is oo 
and f{e,n,k\ ,ky) is an expression in terms of F-functions with arguments in terms of integer 

' However the amount of moments is far too low to reconstruct the general expressions depending on n using the 
method outlined in lElll. 
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linear relation in the parameters n,ki,. . . ,^v-i and e might occur linearly in the form re with r 
being a rational number. 

Then given such a sum, we will present tools to derive the required coefficients in terms of 
indefinite nested product-sum expressions. In particular, given this format, the coefficients can be 
expressed, whenever possible, in terms of harmonic sums, S'-sums and more generally in terms of 
cyclotomic sums [^] defined as follows: 

^(ai,bi,ci),...,(as,b,,cA{^ij ■ ■ • i^s\n) = ^ - ; , , s,,, V 1 ! , , \,., • • • ~, '■ , , \,. (2-2) 

with ai,bi,Cj G N (a,- > bj, Ci 7^ 0) and Xi being an element from the underlying field lEC (containing, 
e.g., the complex numbers); note that choosing a,- = 1, bj = and Xi £ {1,-1} restricts to the class 

of harmonic sums written in the format Sxi^ ,...,x,c, {n) = jr L,-2=i jr ■■■ L/=i ^ ■ For ^ formal 
but lengthy definition of the more general class of indefinite nested product-sums w.r.t. n we refer 



to [32 1 . This class includes indefinite nested sums as in where also hypergeometric terms (like 
binomials, factorials/F-f unctions, Pochhammer symbols) might occur as polynomial expression in 
the numerators and denominators of the nested summands. 



Example 1. Consider, e.g., one of the simple sums arising within the problem of [23]: 



F( ) = V " V ^ ^2^^'^'^(-iV'a2+i)n2-e)r(f+2)r(-f)r(-f+yi+4)r(-ji+n-2)r(e-7i-j2+«-5) 

• n ■ n (e-10)(e-8)(e-2)er(f-e)r(£±5)r(f+„+l)r(-j,-i2+n-4) 

;i=o j2=o s ^ ^ \ ^ ^ ' . 



^ =f{n,h-h) 

= F_3£-3 +/7_2£-2 + /7_i£-i +Foe^ + ...; (2.3) 



note that (2.3) is not indefinite nested w.r.t. n, since the inner sum (at least in this representation) 
cannot be written in the form YJj2=i2f'^{h) where Z2 G N and f2{j2) is free of j\ and n and since in 
addition the outermost sum cannot be written in the form ^7,=/, fiiii) where Zi G N and f\{j\) is 
free ofn. Then loading the Mathematica packages 

in[ii:= << Sigma. m 

Sigma - A summation package by Carsten Schneider © RISC-Linz 

in[2]:= << EvaluateMultiSums.m 

EvaluateMultiSums by Carsten Schneider — (c) RISC-Linz 

(and defining f as the summand of our sum) we can compute the coefficients {F_3,F_2,^-i} with 
in[3]:= EvaluateMultiSum[f, {{j2,0,n ji 6},{ji,0,N-5}},{n},{5},ExpandIn {e, 3, 1}] 

16(-1)'' (Sn^ + 12n + ll) 16 (n^ + 61J - 6n^ - 80n^ - 81n* + USn^ + 274n2 - 4n - 96) 
135(n+l)(n + 2)2(n + 3)2 45(n- 2)(n- l)2n2(n+ l)(n + 2)2 (n + 3)2 

16(n2-n-8) ° 1 8(n2-n-8) ° 1 2(-l)°(l87n+ 127) (3n2 + 12n+ ll) c 

45(n- l)n(n + 2)(n + 3) i' 45(n- l)n(n + 2)(n + 3) Y ^ 2025(n+ l)2(n + 2)2(n + 3)2 

2(l7n^-231n^ + 121n* + 2063n^-1458n2 -2432n + 960) 16(-l)''(3n2 + 12n+ ll) ° 1 
^' 675(n-2)(n-l)2n2(n+l)(n + 2)(n + 3) 135(n+ l)(n + 2)2 (n + 3)2 ^ ^ i 

2(43n" + 112n"+263n"-216n^-11309n''-16476n^+55837n"^+78164n'^-9B178n^-116688n^+51784n^+30624n-23040) 

675(n-2)-i(n-l)%a(n+l)^(n+2)^(n+3)^ J ■ 

The input {n},{5} means that we suppose that n > 5; this specification of the parameter range 
is necessary for the internal calculations. Loading in addition J. Ablinger's Harmonic Sums 
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package /|7]/ would automatically transform the output sums to the harmonic sums Si{n),S2{n). 
The constant term Fo{n) (which would be too large to print here) introduces in addition the sum 
S'i{n) and the zeta-value ^2 defined by the Riemann zeta function Q = YilLiZ^'' the total timings 
(including the constant term) are 4 minutes. 



To get these coefficients, we expand tlie summand of (|2JJ) in the parameter e 

f{e,n,ki, . . . ,ky) = fu{e,n,ki, . . . ,ky)e" + fu+iie,n,ki, . . . ,ky)e"+'^ + . . . 
by using formulas such as r{n +! + £) = '"^"^"'^I^tj^'' with e = re for some r G Q and 

Bin,\ + -e) = lexp ( £ 1 = 1 f ...^i(n) 

" \k=i J f^k=o ^-^^ 

and other well-known transformations for the F-functions. Finally, under the assumption that the 
sums are uniform convergent, in particular if one can exchange summation signs and differentia- 
tion, one gets the coefficients 

Ll(fi) Ly(n,ki,...,ky-i) 

Fi{n)= ••• L fi{e,n,ku...,ky) 

k\=h ky=U. 

of the expansion (|2.1[). Now the hard task is to simplify these sums further. Here we rely on the 



following summation paradigms [30] in the context of difference fields. 



Deriving recurrences: Given an integer d>Q and given a sum 

L(m,n) 

F{m,n) := £ f{m,n,k) (2.4) 
k=i 

where / G N, f{m,n,k) is an expression in terms of indefinite nested product-sums w.r.t. k, and 
where m = (mi , . . . , m^) and n are discrete parameters; find 

aQ{fh^n)F {fh,n) -\ \- aij{m,n)F {m,n + d) = h{fh,n). (2.5) 

where ao{fh,n),. . . ,ad{m,n) are rational functions in m and n, and h{m,n) consists of an expres- 
sion in terms of sums of the type as F{m,n), but with simpler summands (i.e., less summation 
objects or less nested summation quantifiers). 



Remarks. (1) The underlying algorithms [gJJ utilize creative telescoping that has been originally 
introduced with Zeilberger's algorithm [35] for hypergeometric summands. 
(2) The upper bound L{in,n) might depend integer linearly on n,m or might be 00. In this particular 
case, in order to get ( |2.5| ), limit calculations are necessary using asymptotic expansions of the 
summand objects. For the expansion of harmonic sums, S'-sums and cyclotomic sums efficient 
algorithms are developed [||] and implemented in Harmonic Sums [|I]]. 



Now suppose that we are given a recurrence of the type (2.5) where we succeed in representing 



h{fh,n) in terms of indefinite nested product-sums w.r.t. n. Then we can proceed as follows. 



Recurrence solving: Given such a simplified recurrence, find all solutions that are expressible 
in terms of indefinite nested product-sum expressions w.r.t. n. 



Remark. The underlying algorithms, see e.g. [^, ^], find indefinite nested sum expressions that 
are highly nested: the maximal nesting depth will be the recurrence order. Thus the simplification 
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to expression with minimal depth is a crucial step [32]. In addition, the occurring sums are reduced, 



i.e., they are algebraically independent among each other [|33|]. 



Combination of solutions: Given these simplified solutions, combine them such that the evalu- 
ation agrees with the original sum F{in,n) for the first, say n = 1, + d — I initial values. 



If this is possible, it follows -up to some mild side conditions- that this combination of indefinite 
nested product-sums is equal^ to the definite input sum for all n>l. 

Example 2. The goal is to find a representation in terms of indefinite nested product-sums w.r.t. n 

for the coefficient =/(«j'i) 

, ^ ^ 

n-5 n-ji-6 



^ 135(n+l)! Lj ,• _|_ ,• _ „ _|_ 5 ^ " ^ 

;i=0 ;2=0 ilT^i2 " T J 



in the expansion of sum (23). In the first round, we transform the summand /(n, j'l) to indefinite 
nested product— sums w.r.t. ji with the algorithms implemented in Sigma. First, we compute a 
recurrence for f{n,j\ ): 

i-ji +n- 5)f{n, ji ) + {ji-n + 4)f{n, ji + I) = ji - n + 5. 
Next, we solve the recurrence and get the full solution space L = {c{—ji + « — 4) + (—71 + n — 
4-)Y,jLi -3+n-i ^ ^J" which means that f{n,ji) € L. Now we have to combine the solutions (i.e., 
to determine c) so that they agree with f{n,ji ) at, e.g., ji = 0, i.e., /(«, 0) = Ly2=o ^ ;2^-n+5 • ^S'^^^' 
we turn this sum to an indefinite nested sum (by computing a recurrence, solving the recurrence and 

combining the solutions) and obtain f{n,0) = ''(n-'i){n-2)(n-\)n — ^ + (4 — «) 7- Choosing 

c = f(^n,0) /{n — A) finally yields 

f{n,n) = {-n+n-A)Y^ +^ („_3)(„_2)(n-i)„ + (ji-n + 4)X^. 

,= 1 -> -r" t I 

With this indefinite nested sum representation w.rt. ji of the summand f{n,ji), we start the final 

round: we compute a recurrence for ([2.6|), solve the recurrence and combine the solutions to get 

the coefficient F^2{n) computed in Example 0. 

In a nutshell, given a multi-sum of the form 

L{fn,n) L(ih,n,k) Ly{m,n,k,k[ i ) 

F{in,n)= XI L ••• L f{m,n,k,ki,...,ky) 

k—l kl—l\ ky—ly 



f(ih,n,k) 

where / itself is an expression in terms of indefinite nested products-sums (in particular the class 



of sums given in (2.1) is covered), we apply the following method to transform F{fh,n) to an 
expression in terms of indefinite nested product-sums. 

1. Transform the outermost summand f{m,n,k) to an expression in terms of indefinite nested 
product-sums w.r.t. k by applying the method recursively to all the arisings definite sums 
(i.e., the parameter vector in is replaced by («,m) and the the role of n is k). Note that the 
occurring sums in / are simpler than F(n) (one definite sum less). If the summand is free of 
sums, nothing has to be done. 

^We emphasize that in each computation step of this transformation we produce proof certificates which enables 
one to verify rigorously the correctness of the result. 
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2. Compute a recurrence (2.5) for the sum (2.4); note that by construction the right hand side 
might be again an expression in terms of definite muhi-sums, but these summands are of 
simpler format than /. Apply the method recursively to these sums and compute a right hand 
side representation which consists of indefinite nested product-sums w.r.t. n. 



3. Solve the recurrence (2.5) in terms of indefinite nested product-sums w.r.t. n. 



Compute d initial values (i.e., specialize the parameter n to appropriate values from N, say 
n = l,l + l,...,l-'rd — I, and apply the method recursively to the arising sums where m\ 
takes over the role of n and the parameters are in. If no parameter is left, the expression 
is a constant. In particular, if there is no sum left, nothing has to be done. Otherwise, in 
the outermost summations the upper bound is oo. In this case, our method is applied to 
transform the summands to indefinite nested product-sums w.r.t. the outermost summation 



index. Usually, the derived sums can be transformed to multiple zeta values [19] or infinite 



versions of 5-sums or cyclotomic sums [Q]. Otherwise, they are kept in the given format. 

5. Try to combine the solutions to find an indefinite nested product-sum representation w.r.t. n 
oiF{m,n). If this fails, ABORT. 



Remarks. (1) The existence of a recurrence in step 2 is guaranteed by using arguments form [27] 
and [^5|]. Only computation issues are a bottleneck. Usually, we succeed in finding recurrences 
when / consists up to 100 indefinite nested product-sum objects. If / is more complicated (or if it 
seems appropriate), the sum is split into several parts and the method is applied separately. 

(2) Termination: The method is applied recursively to sums which are always simpler than the 
original sum (less summation quantifiers, less parameters, or less objects in the summand). After 
finitely many recursion steps, one arrives at the base cases where no summation quantifiers arise. 

(3) Success: If the method does not abort in step 5, it terminates and outputs an indefinite nested 
product-sum expression w.r.t. n. As a consequence, finding not sufficiently many solutions of a 
given recurrence of the type (E3) in step 5 is the only reason why the method might fail. For general 



multi-sums this failure would happen all over. However in the context of Feynman integrals, we 



found almost always d linearly independent solutions of the homogeneous version ( |2.5| ) and one 



particular solution of the recurrence itself; in these cases, the solution space of ( |2.5| ) is completely 
determined, and the failure in step 5 cannot occur. The (for us surprisingly) rare case that not 
sufficiently many solutions are found was always an indication that the sum representation of the 
Feynman sums could be improved so that afterwards our method worked. 

(4) Subtle details: The input sums coming from Feynman integrals are rather challenging. In 
each step of the method, pole issues arose (lower and upper bounds must be updated during the 
calculations and thus compensating terms must be computed separately, the initial values must be 
chosen carefully, etc.). As a consequence, the 5 line method above implemented in the package 
EvaluateMultiSums requires currently about 8000 lines of code in Mathematica (not counting 
the implementations for recurrence finding and solving which are part of the package Sigma). 

The package EvaluateMultiSums has been applied successfully to various complicated 
3-loop ladder graphs [^, ||]. In particular, the graphs from [Q] could be calculated by reducing it to 
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several multi-sums, one of them being 

"v V V (a!)(^)0;l)(-^';''-^)(-'n-^-^)("""-r'-^)'-H-'+"-'/-^-.v-3)!(^-i)! 

2-1 2-1 2-1 2-1 2-1 2-1 {-l+n-q-2)\{-i+n-\)(n-q-r-s-2)(q+s+\) 

j=0k=0l=0 q=0 5=1 r=0 

(- 1 ) -i+^-'+«-'?-3 \ASi{-j + n-\)~AS,(-i + n~2)-2S,{k)~{S,{-l + n~q~2) 
+ Si{-l + n-q-r-s-'i)-lSi{r + s))+lSi{s-l)-lSi{r + s) 

3. Mass production 

For various problems, like e.g. [^, ^, several thousand sums have to be evaluated. In particular, 
for the brand new calculation of fermionic contributions to the gluonic massive operator matrix 



elements [23] the £-expansion of a 2 GByte expression consisting of 2419 multi-sums was calcu- 



lated; one of the simple sums is, e.g., (|2.3|). With a lot of computer resources this problem could be 
tackled ^. However, we can do it much better by using the following routines which are available 
in the new package 

in[4]:= << SumProduction.m 

SumProduction - A summation package by Carsten Schneider © RISC-Linz 

1. Reduction to key sums: First, we reduce the 2 GByte expression (stored in expr) to key sums 
with the function call 

in[5]:= compactExpr = ReduceMultiSuins[expr, {n}, {5}]; 

The reduced expression compactExpr is only 7.6 MByte large and contains only 29 sums and 
15 terms free of sums; in total it took us 6 hours and 53 minutes to obtain this reduction. 

Remark. Internally, the 2419 sums have been synchronized w.r.t. the occurring summation ranges 
(taking for each class the maximum of the lower bounds and the minimum of the upper bounds). 
As result, we obtained only 4 sums with equalized summation ranges 

n— 5 (2 '1—5 '1—12—5 n— 5 °° 

Y^Y^hi{£,n,i2M), ^ h2{e,n,i2,h), ^ /z3(e,?i,/i), ^ /j4(£,m,/i) 

'2=5'i=0 12=0 11=0 (1=5 11=0 

plus a large term free of summation quantifiers. Next, all the occuring Pochhammer symbols, 
factorials/T-functions, and binomials are written in a basis of algebraically independent objects 



plus the extra object (—1)" (if necessary); for details see [^, Sec. 9] and [29, Sec. 6]. Finally, 
the expressions are split further to get the form Y.h{n, (/2,)/i,£) *r{n, (/2,)/i,£) or h{n,e) *r{n,e) 
where h stands for a (proper) hypergeometric term in n (and ii , 22), i.e., being a product of binomi- 
als/factorials/Pochhammers in the numerator and denominator, and r{n, (/2,/i),e) being a rational 
function in n, e (and 11,12)', note that r might fill several pages. 

2. Computing £-expansions in parallel: Next, we compute the coefficients of the £-expansion as 



outlined for sum (|2.3|); as it turns out, the time to calculate the expansion of a key sum is similar to 
calculating just one typical candidate within the 2419 sums contributing to the corresponding key 
sum. In order to produce these expansions automatically, we developed the routine 



For testing we calculated around 800 of these sums with the help of 16 Mathematica processes in around 4 days. 
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in[6]:= ProcessEachSum[compactExpr, {n}, {6},ExpandIn {ep, 3,0}] 

which sequentially applies EvaluateMultiSum with the corresponding input parameters to the 
occurring multi-sums in compactExpr. This step took in total 2 hours and 35 minutes. 

Remark. Internally, it takes the first multi-sum and generates a file with the name SUMl. If the 
result is computed, the file is updated with the result. Then the routine continues with the second 
sum provided the file SUM2 is not existent on the hard disk. In this way, ProcessEachSum can be 
executed in parallel for mass productions (for even larger problems than this). 

3. Combining the subresults: Finally, the coefficients of the expansions of the subresults are read 
from the hard disk and are summed up to the final coefficients of the expansion: 

in[7]:= result = CombineExpression [compactExpr, {n}, {6}]; 

Remark. Internally, the expressions are reduced further by eliminating all algebraic relations of the 
occurring sums and products [[T^, ^ |3|]. This step took 21 seconds. 

The final result can be expressed in terms of 

C2 , ?3 , (- 1 )" , -^l (") , -^2 («) , ^3 («) , ^2, 1 («) , ^3 J («) , 52,1, 1 («) 

and requires about 100 KByte memory. The total calculation took around 9 hours and 30 minutes. 

In conclusion, we are ready to go into mass production for various challenging problems being 
of similar type as outlined in [^, ^, ^, ^ ]. 

Acknowledgment. This work has been supported in part by DFG Sonderforschungsbereich Tran- 
sregio 9, Computergestiitzte Theoretische Teilchenphysik, Austrian Science Fund (FWF) grant 
P203477-N18, and EU Network LHCPHENOnet PITN-GA-2010-264564. 
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